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1* Introduction 

In this paper we describe procedures for computing linear and nonlinear aspects of acoustic wave 
propagation in the atmosphere. Acoustic waves originate from a source that has compact support in 
space, as well as in time. Such a phenomena occurs in a situation such as shuttle take off or any general 
explosion. The current work concentrates on a situation in which the atmosphere is isothermal. 
Although this assumption somewhat over-simplifies the physics, the mathematical and computational 
aspects of the problem remain difficult. The main difficulty arises from the effect of gravity, which is 
included in this model and plays a crucial role in the nature of the propagation. In the acoustics litera- 
ture, similar models were used to investigate explosive waves. We cite references [1] and [2] where 
linear models were considered for such purposes. Our model is an improvement over the model pro- 
posed in [1], Our work appearing in [3] and [4] discuss analysis of the problem in two dimensions. Our 
work appearing in [5] discusses the first computational effort in the three-dimensional situation but only 
the linear case was discussed. The work reported in this paper includes numerical aspects of nonlineari- 
ties in the three-dimensional situation. 

The main goal of this paper is to capture the effects of these atmospheric acoustic waves numeri- 
cally. While the goal seems to be simple, there are several obstacles in the process of solving the prob- 
lem under consideration. The nonlinear problem for the acoustic variables cannot be directly solved 
with a reasonable computational effort unless the problem is somewhat simplified by some techniques. 

In the sequence of work reported in [3] and [4] the nonlinear problem was reduced to a sequence 
of linear problems using appropriate asymptotic expansions. Then one could show the linearized prob- 
lems were well-posed. There was a question that we could not resolve about the asymptotic expansion. 
The expansion may not be uniformly valid. Here, we provide partial justification by considering an 
appropriate one dimensional model. We compute the linearized solutions and show the sum of the first 
two terms of the expansion yield indeed the solution of the nonlinear problem. Unfortunately, this can 
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be done only in a one-dimensional situation, as the cost of computing solutions of the nonlinear prob- 
lem directly is extremely high. This is because the actual nonlinear acoustic field contains the mean flow 
fluid quantities which are orders of magnitude higher than the acoustic quantities. In particular, the 
truncation error resulting from the numerical scheme will dominate the acoustic effect unless a very 
fine grid and extremely small time steps are used. In the calculations reported on the one-dimensional 
situation, 400 mesh points and about 10,000 time steps are needed to directly compute the nonlinear 
problem. To carry out this procedure to higher dimensions, it is clear that one requires large computa- 
tional requirements. The justification for the validity of the asymptotic expansion that we provide here 
is computational, and a mathematical proof is a hard task as pointed out in [3]; therefore we leave it as 
an open question. 

The plan of the paper is as follows. First, we shall present the governing equations and present 
the problems that we intend to solve. In these problems, boundary conditions play a crucial role. The 
boundary conditions that we use here are different from and better than the ones we used in [5]. The 
derivations are not straightforward. They are presented in the appropriate sections. We devote section 2 
of the paper to the discussion of the one-dimensional model where the problem can be reduced to a 
scalar problem. In section 3, we discuss the sequence of linearized problems and in section 4, we dis- 
cuss the nonlinear system in the one-dimensional situation. In section 5, we discuss the axisymmetric 
three-dimensional problem. We also show the numerical results in sections 3,4 and 5 arising from our 
theory and computations. In section 6, we show the numerical schemes that we used for our computa- 
tions, along with the numerical treatment of our boundary conditions. 


1.1 Governing Equations 


We shall begin with the statement of the fluid flow problem that governs the acoustic phenomena. 
If P* is the ambient pressure and h is the scale height, then the non-dimensional form of the Euler 
equations (the equations of continuity, balance of momentum, and energy) are 


jfc + d/v(pq) = 0, 


-|a + (q. V)q = -J-Vp- ±k, 
at YP Y 



l£- = 


P y 1 


( 1 . 1 ) 

( 1 . 2 ) 

(1.3) 


_ 1 ~ 

Note that, in equation (1.2), the forcing term k (k is the unit vector in the z direction) arises due 

Y 

to the forcing term per unit mass -gk in the original variables which is due to gravity. In equation 
(1.3), f (x ,y,z,/) dictates the space-time dependency of the source, and e measures the energy release 
per unit volume. For the case of an instantaneous energy release, e is given by 



where Q 0 is the total energy released at time t = 0. The initial conditions are 


p = p = e 1 , q = 0 at t = 0, 

(1.5) 

which represent a calm atmosphere and an exponentially decaying pressure and density. 

The boundary condition at z = 0 is 


■P 

ii 

o 

(1.6) 

which states that the vertical component of the flow is zero at z = 0. 


The acoustic expansion is based on e << 1 and represents the flow as small changes superim- 
posed on the flow of the ambient state. We note that the ambient velocity is zero, but pressure and 
density have the form e ~ z . Thus, the expansions are 

q = eu + e 2 ui + • • • , 

(1.7) 

p = e~ z + ep + e 2 /?! + • • • j 

(1.8) 

p = e~ l [l + ea + e 2 ai + • ■ • j 

(1.9) 

where u = («,w) and uj = Quantities w, u\ and w, Wj are the r and z components respec- 

tively of the acoustic velocities, in cylindrical coordinates. We substitute expansions (1.7) -(1.9) into 
equations ( l.l)-( 1.3) , initial conditions (1.5), and boundary conditions (1.6) and retain terms of order 
e to obtain the field equations. Therefore, the field equations are 

u _ 

a, + u r + — + w z - w = 0, 

(1.10) 

1 

«( + —p r = 0, 

(1.11) 

w, + 1 ft - - 0, 

Y Y 

(1.12) 

p, + y u r + yw z + y— - w = f (r ,z ,t). 

(1.13) 


where y = 1.4 and /(r,z, t) describes the nature of the acoustic source and r = V x 2 +jT 2 . 

Boundary Conditions are: 

w(r,0,0 =0, t> 0 (1.14) 


u(r,z,/)/r is bounded as r— >0. 


(U5) 
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Initial Conditions are: 


p = a= u = w = 0 for t - 0. 


(U6) 


In addition, we also need a radiation condition for the problem that describes the behavior of the wave 
at infinity. In this problem, a radiation condition must be included for the uniqueness of the solution. 
The condition here is that the waves decay at infinity and, when imposed at finite distances, take the 
form 

(u -p) t = w z -w + — at r = R , (1.17a) 

y Y r 

(w - ^-p) t = ■!■(/> - a) - - — - u r ) at z = L. (1.17b) 

Y Y y r 

The derivation is based on [8] and is given in section 5. 

We shall denote the problem formed by equations (1.10) - (1.17) by (P). Our procedure advo- 
cated in [3], as well as in [4], is to solve for the next correction and add the correction to the solution 
of (P) to obtain the solution of the nonlinear problem. In this process, the equations governed by the 
next order terms contain exactly the same differential operator, except the right hand side contains the 
solution of problem (P). We call the the resulting problem of the next order to be (P’) and is as fol- 
lows: 



(U8) 

(1.19) 

( 1 . 20 ) 
( 1 . 21 ) 
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The boundary and initial conditions have the same form as before and are given below. 

Boundary Conditions are: 

wjfr.O.f) = 0, t > 0, (1.22) 

u x (r,z,t)lr is bounded as r-> 0. (1-23) 

Initial Conditions are: 

p l = (7 1 = = W] = 0 for t = 0. (1.24) 

The radiation condition takes the form 

/lx / 1 1 

(“i - —Pi)t = -(-wi, + — wj - — ) + —ap r - uu r 

- wu z + — [p(w-y— ) - Y p(u r + w 2 ) - ( up r + wp z )] at r = R , (1.25a) 

Y r 

( W 1 “ = “~( w l - y— - Y«1 r) + -“(Pl - ^l) + ^ +W * + ~ CT (P - 

y Y r y y Y 

- —P(w “ Y“) + P( u r + + — (w/? r + wp z ) at z = L. (1.25b) 

y r Y 

It is shown in [3] and [5], these linearized problems are well-posed. However, it remains to be deter- 
mined if the expansions proposed in equations (1.6) - (1.8) are uniformly valid. If they are, then we 
have a justification to assume that the full nonlinear problem is well-posed in the sense of the results 
reported in [3] and [5]. To demonstrate the validity of the asymptotic expansion, we consider a 
sequence of one-dimensional models. We show that the resulting linear problems from the asymptotic 
expansions indeed lead to the solution of the original nonlinear problem. At the outset, let us 
emphasize that even the linear one-dimensional problem does not have a closed form solution. For this 
reason, we solve the problem in two essentially different ways to reach an agreement in the solution 
that we obtain through numerical methods. In the next section, we reduce the problem to one that is 
governed by a scalar wave equation. While the field equations become simple, the formulation of the 
radiation condition becomes a laborious task. The derivation is included. In section three, we consider 
the linearized systems. In section four, we solve the nonlinear equation. Numerical methods and 
derivation of the boundary conditions are included. Here we compare the results obtained in the 
preceding two sections. In section five, we solve the three-dimensional axisymmetric model. Again, the 
boundary conditions are paid special attention. Also, we formulate the problem in the absence of grav- 
ity and show the derivation of appropriate nonreflecting conditions. We solve this problem and compare 
the effect of gravitation with the solution obtained from the original systems. In section six, we con- 
clude the paper with the numerical methods that we used. 
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2. One Dimensional Model Problem for the Sound Propagation in the Atmosphere. 

In this section, we propose a method of solving the one-dimensional model of sound propagation 
in the atmosphere. As mentioned earlier, the radiation boundary condition that is suitable for the 
numerical calculations is a difficult task. We derive the boundary conditions. The propagation is res- 
tricted in the vertical direction where the effect of gravity plays a crucial role. For computational pur- 
poses, the distance which extends to infinity needs to be truncated. Let us denote this finite interval to 
be [0,L ]. The ideas behind deriving the radiation condition here are very similar to the one proposed 
recently in [6] in which wave propagation in certain types of viscoelastic materials are considered. The 
key idea here, as in [6], is to construct an interface problem with an interface at z = L , where, typi- 
cally, a well suited radiation condition that is best for numerical calculations can be posed. For z > L , 
the problem will be solved by an asymptotic procedure that is valid for small times. This procedure is 
similar to the one reported in [7]. But, it turns out that these boundary conditions are not necessarily 
accurate only for smaller time scales. To present the idea, let us introduce the governing equations 
obtained from equation (1.6)- (1.8). The equations are 


a t + w z - w = 0, 

(2.1) 

w < + J p i - J (P-®)= 0. 

(2.2) 

p, + yw z - w = f. 

(2.3) 

These equations are defined over 0 < z < . Also, these equations, unlike in the higher dimensions, 

can be reduced to a single scalar equation as follows 


+ w, = Y (/ - f z ) = G(z,r). 

(2.4) 

Boundary conditions are 


p 

ii 

p 

(2.6) 

Lim w(z,/) is bounded. 

z —-><*> 

(2.7) 

Initial conditions are 


w(z,0) = 0, 

(2.8) 

w,(z,0) = 0. 

(2.9) 


Equation (2.4) can be simplified by the transformation 


w(z,r) = e Xz v(z,r), 


with a choice of X = as follows: 



( 2 . 10 ) 
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1 “7 Z 

V« - V„ + —V = <? G(z,0 = g(z, 0- 

Suppose we want to impose boundary condition, exact or otherwise, at the finite distance L that 
simulates the behavior at infinity. We claim that these boundary conditions are of the form 

v z = Bv + Kg , (2.11) 

where 5 and are operators to be determined. In the event that g has a compact support in 
0 < z < L , the boundary condition will take the form 

v z = Bv. (2.12) 


2.1 Procedure 

Initial conditions suggest that we consider the problem in the Laplace transform domain. Doing 
so, we obtain 

s 2 v = v zz - JV + g. (2.13) 

Here v = v(z,j) and g = g(z„y). The main theme here is to consider an interface problem as fol- 
lows 

V„ = V zz - JV + g(z,t), 0< z< L , (2.14) 

which is the given field equation. Now define U (z,.s;r) such that 

s 2 U = U zz - jU + g(z,t), z> L, (2.15) 

with 

U (L ,s;t) = v(L ,/)• (2.16) 

We remark that here t is considered as a parameter. Now it is easily verified that 

v(z,s) = J e~ s ‘U (z,s;t)dt. (2.17) 

0 

The plan here is to solve equation (2.15) subject to the interface condition (2.16) and a radiation 
condition which insures the boundedness of v(z,r) , which translates in the transform domain to 
requiring the boundedness of U{z,s\t). As mentioned earlier, our procedure for constructing the solu- 
tion will be an asymptotic one with s as a parameter and for large s , reflecting the fact that solutions 

will be accurate for smaller times. However, asymptotic methods are known for their power. Even 

though the expansion is valid for small time in the theoretical treatment, need not be small for real 
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applications. To elaborate this statement, one may consider an asymptotic procedure given in [7], which 
in theory was valid for a parameter a approaching infinity, but it turned out surprisingly enough to be 
valid for a = V 5 in numerical experiments. Thus, the results are anticipated to be accurate for more 
than just small times. 

To solve equation (2.15), we shall use elementary concepts, except incorporating appropriate 
asymptotic ideas. We first seek a particular solution. Let W (x >s;t) satisfy 

Wzz - ( S 2 + J)W = -g(z,t). (2.18) 

Then W may be sought in the form 

~ Wniz) 

W(z 9 s;t) = S — V" - (2.19) 

k = i * 


Substituting the expansion in equation (2.18), we obtain 


) _ s2 (— r + ~r + ~~r + 


. iw 2 w 4 

} 4 ( ^ + ^4 + 


Equating powers, we obtain the following: 

{power 5°) W 2 = g, 

(power s~ 2 ) W 2 "~ W 4 - jW 2 = 0, W 4 = g" - jg, 

(power s~ A ) W A "-W 6 -\w A = 0, W 6 = g ^ - \g” + ~g, 

4 L lo 


etc.. 

Remark: 

One may consider odd powers in the above expansion and verify that all corresponding terms are 

zero. 

Thus, 



Let w = U - W (i.e., the homogeneous solution ). Then w satisfies 


- (s 2 + — ) w = 0. 


To solve the homogeneous equation, we make some intelligent choices. We incorporate the fact that as 
z approaches infinity the homogeneous solution should vanish. In particular it should be outgoing. This 
lead to the expansion 
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w 


e -s(z-L)( Uo+ El + 
s 



( 2 . 21 ) 


Substituting the above expansion in equation (2.20), we have 


-mu.-* of * v 4 ♦ * *£ * ...) - >. ♦ ^ ^ ♦ ...) = o 


u 


Ui 


Ux 

s 


U 7 


Recall that w = U -W. We also have the requirement that U(L,s;t ) = v(L,t). Thus, 
w(L,s;t) = v(L,t)~ W(L,s;t), indicating that the first term independent of s on the right is 
v(L ,t). Using this fact in the previous expansion (2.21) and by equating the s 1 powers, we obtain 

' U 0 ' = 0; U 0 (L,t ) = constant-, U 0 (L,t) = v(L,r). 


The next term we obtain in the expansion is 

-2 US + U 0 " - jU 0 = 0. 

Again we note that at z = L the right hand side of w = U — W does not have inverse power of 
s , yielding U X (L ,t) = 0. Thus, integrating for U i, we get 

Ui(z;t) = ~jv(L,t) (z-L) 

The next term in the expansion has some interesting features. This time, forcing terms will appear in 
the expansion. When doing the expansion, we find that 

— 2U 2 + U x " - jU x = 0. 

We see as before from the construction of w that 


U 2 (L,t) = -W 2 = -g(z, /). 

Thus, integration yields 

U 2 = -^-v(L,t)(z-L ) 2 - g 


A similar procedure yields 


U , = 


128 


v(L,t)(z-L ) - 


3182 


v(L, 0 (z-L ) 3 - g'Xz-L). 


We note here that on the boundary z = L , C /3 = 0. Moreover, the remaining terms in the 
expansion, when evaluated at z = L , yield only terms involving g and its derivatives. Thus summar- 
izing the results, we have 


U (z,s;t) 


JL 
« 2 


+ e 


-^- L) (U 0 + — + ^- + 


) + 


( 2 . 22 ) 



- 10 - 


Differentiating with respect to z, we have 

a £ z , t \ U lt U 2 z / i \ C/ 1 U 2 

U z (z ,s\t) = 4 + ^- j(z “ L) (C/o,z + — ^ + -T- + ...) - ^ ( "“ L) (t/o+ — + -r + ...). 
s l s s l s S L 


Using the relations that we obtained for U Qt U U 2 , we obtain 

U z (L,s;t) = -sv(L,t)- l v(L - ^ + M 1 ’*) +F(g,s ), 

O 5 1 J 


(2.23) 


where f is a functional of g and involves the inverses of s 4 and higher powers. Neglecting F and using 
the definition (2.17), we have 


v,(L, s) = -SHL. S ) - + 

OS s 


(2.24) 


Inverting the Laplace transform, we obtain our desired boundary condition as follows: 


i(L .0 = -v,(L.O - iJv(L,t)rfx + Jg(L,T)dT, 


(2.25) 


t (L ,t)+w,(L ,t) = - -t-J w(L ,x)dx + J s(L ,x)e 2 dx. 


(2.26) 


Thus for computational purposes, we shall replace the boundary condition (2.7) by the approximate 
boundary condition (2.26). 


3* Verification of the Asymptotic Theory and the Numerical Procedure 

As pointed out in the first section of the paper, it is difficult to construct an analytic solution, 
even for the linear problem. Thus we content ourselves with numerical solutions. To verify the correct- 
ness of the linear problem, we compute the results through two different formulations. One is to solve 
the reduced wave equation (2.6) subject to the initial conditions (2.8) and the boundary conditions 
(2.6) at z = 0 and approximate radiative condition at (2.26) at z = L. Let us emphasize here that we are 
solving the problem over the interval [0 ,L ] and not over [0,°°). For the ease of the reader, let us 
restate the problem consisting of the reduced wave equation which is indicated below and we shall call 


it (P-W) for future reference. 

w,t - w zz + w z = — (/ - f z ) - G (z ,/). (3.1) 

T 

w(z,0) = 0, (3.2) 

w.(z,0) = 0, (3.3) 

w(0,O = 0, (3.4) 

1 1 ' ' L 

w z (L ,t)+w,(L ,t) = —w(L ,i) - — f w(L ,x)dx + jg(L,T)e 2 dx. (3.5) 

L 8 0 0 
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Solving (P-W) numerically is straight-forward. It is solved by an explicit finite difference formula. 
This is a standard discretization. The integrals in the boundary condition (3.5) are evaluated by a rec- 
tangular quadrature formula. It could be shown that problem (P-W) is well-posed in the sense of a 
bounded growth of energy. Related ideas can be found in MacCamy [6]. However, the situation here is 
different and the related theoretical results will be reported elsewhere. 

As mentioned earlier, the comparison could be made by solving the system (2.1) - (2.3) with 
appropriate initial and boundary conditions as described in the last section. Still the problem is to be 
considered in the truncated domain [0,L]. However, the system involves the primitive variables for the 
acoustic components p , p and w. The boundary condition (2.26) thus need to be recast in these vari- 
ables in an appropriate manner which is suitable for numerical computations. Again, we summarize the 
equations that we wish to solve. 


d, + w t - 

w = 0 

(3.6) 

+ 

*< 1- 
;? 
i 


(3.7) 

Pi + Y w z ~ 

w = f 

(3.8) 


The initial conditions are 

p = a = w = 0 (at t ~ 0). (3.9) 


The boundary conditions are 


w(0,0 = 0 (3.10) 

and (3.5). We replace w z in (3.5) using equation (3.8) so that the left side of the boundary operator 
becomes the rate of change of incoming Riemann variable. The result is 


— (/> - yw) = 


(1 


~ ') w + / + •g'J W ( L “ Yj 8( l ,i)e 2 dx. 


(3.11) 


We shall call the problem (3.6) - (3.11) by (P-Ll). We briefly show the calculations that identify the 
Riemann variables associated with the system. First, we note that that the system (3.6) - (3.8) can be 
written in the form 

u, + Au* = H (3.12) 


with 


u r = ( a,w,p ), 

H r = (H'.ytp- a),w + f), 


(3.13) 


(3.14) 
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and 


A = 


0 1 
0 0 
0 y 


0 

y 

o 


(3.15) 


It is easily shown that the eigenvalues of A are (1,-1, 0). Corresponding left eigenvectors are 

ei = (0,Y,l) r , 
e 2 = (0,-Y,l) r , 
e 3 = (- 7 , 0 , 1) T . 


Multiplying on the left of (3.12) by each of these vectors, we obtain 


R t + R z = w + p - a + / , 

(3.16) 

S t - S z = w-p + a + /. 

(3.17) 

T t = w - yw + /. 

(3.18) 


Here 


and 


R = p + yw 


S = p - yw t 


are the outgoing and incoming Riemann variables respectively. T = p - ys carries no information, as 
it has zero speed. These variables play a crucial role in the numerical implementation of boundary con- 
ditions. For example, when we implement the boundary condition (3.11), £ is determined at the 
current time level. However, to obtain the value of p and w\ we use equation (3.16) (discretizations 
may be seen in section 6.). This determines R . Once R and S are determined, p is determined by 

P = (2 = I), (3.19) 

W = (2 = L). (3.20) 

A similar consideration is given at z = 0, where w = 0. To determine /?, one can use equation 
(3.17) for S . Once S is determined at z =0, one immediately obtains 

p = S (z = 0). (3.21) 

Figures 1. and 2. show a comparison of the results obtained through solving the problems (P-W) 
and (P-Ll) numerically. For the test case presented here, the source function / was chosen to have a 
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compact support in space and time as follows: 


/(*. 0 


-a(z - zq ) 2 . 

e sin 2 tc t 


0 < t < T , lz - z 0 l < e 


otherwise. 


(3.22) 


Here a = 40, z 0 = 0.25L , e = 0.025L and T = 1. The results are the same. The results shown 
here are after 1000 and 7000 time steps, respectively, and the pressure distributions are shown in Fig- 
ures 1. and 2. In Figure 2., there is a little difference which is due the boundary conditions. However, 
the scale is small and the difference is below the truncation error. Thus, with this agreement, we con- 
tend that the results are nearly equal to the analytic solution of the problem. In Figure 3, we show the 
behavior of the pressure wave as a function of both space and time and time reaching up to 10,000 
steps or 10 periods. We see that the wave moves in two directions: one towards the radiation boundary 
and the other towards the ground. As one should expect, the wave that hits the ground reflects at 
z = 0, and, eventually, all decay to zero. In the calculations reported here, A z = 0.025, L = 10, 
and the perturbation parameter e = 0.005. There were 400 grid points used for the computations, indi- 
cating a fine grid. Thus, the results that we obtained here had to be thought of as exact. 


3.1 Nonlinear Correction and the Fully Nonlinear Solution 

Now we proceed to investigate the nonlinear effects. First, we propose a second-order correction 
to the linear problem. This problem results by retaining 0(t 2 ) terms as in equations (1.18) - (1.21). 
These equations are 

a u + w u - Wj = -( Gw) z + aw = f u (3.23) 

1 (P i _CT i) 

wi f + —Pn = - [aw, + ww z ] = / 2 , (3.24) 

Pu + ywiz - W X = Gf + GW - ypw z + 7 ww t = / 3 . (3.25) 

Initial conditions are 

Pi = G\ = Wi = 0 for t = 0. (3.26) 

The boundary condition at z = 0 is 

w ! ( 0, r ) = 0, t > 0. (3.27) 

The radiation boundary condition at z = L has the same form as before, which is 

t t \_ L 

fi+ ■£[ w^L^d-t - y\ e 2 H(L,i)d%, 

° 0 0 


jj(p i - ywi) = (i - £) wi + 


(3.28) 
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where 

-\ L i 

H(L,x) = [( p - a) w zz - 2ww zl - (7 + 1 )w z w t ]e 2 - —fc z . 

Let us denote the problem described by the equations (3.23) through (3.28) by (P-L2). This 
problem has the same form as the problem (P-Ll). The same numerical scheme that applies to (P-Ll) 
holds here, including the implementation of the radiation condition. 

Figure 4 illustrates the same situation that was considered in Figure 3. The results shown in Fig- 
ure 4. contains the nonlinear correction. In particular, the nonlinear field is computed according to 

p = zp + e 2 p h (3.29) 

The results given here look roughly the same as the linear case. However, for the linear case, the pres- 
sure profile is smooth. Moreover, the amplitude of pressure differs significantly in the nonlinear correc- 
tion. To validate this nonlinear solution, we solved the fully nonlinear problem in an independent 
manner as described in the next section. Indeed, we obtained the nonlinear behavior correctly by com- 
bining the two linear solutions obtained from (P-Ll) and (P-L2). 


4. Fully Nonlinear One-Dimensional Problem 

To solve the full nonlinear acoustic equations, first we have to obtain the appropriate governing 
equations. It is also important to obtain the appropriate radiation condition for the nonlinear system. 
Let us begin the discussion with the governing equations. 

The field equations governed by the fluid flow are 

Pi + P<?2 + qPz = 0, (4.1) 

<h + Wz + -ZZPz = ( 4 - 2 ) 

YP Y 

Pi + Wz + YMz = e/(z,f)P- (4.3) 

When written in the matrix form, we have 

U, + AU 2 = H, (4.4) 


where 


and 


U = (p ,q,p) T , 

H = (0,— ,e/p) r 
Y 


A 


<7 P 
0 q 
0 yp 


0 

1 

YP 

a 
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The eigenvalues of A are X = q + A/ , q - A/ and q. Left eigenvectors corresponding to 

\ p \ p 

these eigenvalues are 


e, = (0,YVpp,l) r , 
e 2 = (0,-YVpp,l) T , 

e 3 = d,0,--^-) r , 

IP 

respectively. Multiplying equation (4.4) by ei , e 2 and e 3 we respectfully obtain the following equations: 

1 


(Y^PP.l) 

[p] + + V f XtKpp.D (p 

1 = (Y^PP.D 

r a 

, (4.5a) 

(-yVpp.l) 

[?1 + («? - V ^)(-7^pp,i) (; 

l\ = (-Y^PP.l) 

Y 

,(4.5b) 



CO 

l 

If 


(4.6) 


For sources that have a compact support as described in the numerical experiments earlier, equation 
(4.3) suggests along a streamline that 


P =^P Y > 


(4.7) 


where 


A = 


P o 
Po Y 


With the above equations we can describe the fully nonlinear acoustic field and the associated radiation 
conditions. Let 

p = e~ z {\ + p'), 
p = e~ z (l + p'), 
q = u'. 

Thus, rewriting equation (4.4) in terms of the acoustic variables and dropping the "prime" notations, we 
obtain the field equations in the matrix form as follows: 


P 

u 

P> 


0 

1 


u (1 + p) 

0 u 

Y ( 1 + P) 

0 yO + p) u 


r 

P 

u 


M ( 1 + P) 

P ~ P 
Y(1 + P) 

M (1 + p) + e/(l + p) 


(4.8) 
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These equations are to be solved subject to the initial conditions 

p = p = u = 0 at t = 0, (4.9) 

and boundary conditions 

u = 0 at t = 0, (4.10) 

and a radiation condition which we shall derive next at z = L. 

When the above mentioned acoustic expansions are used in equations (4.5a) and (4.5b) together 
with the isentropic relation (4.7), we obtain the following equations: 

1= L ii 1=1 -i=i 

R, + (u+ (1 +p) * )R, = ~j-+ Mu + (1 +p) * ](1 +p) 2Y + SL(1 +p) * (4.11) 

y y y 

S, + (« - (l+p)^” )S, = -±+ Mu - (l+p) V Kl+p)^ - -Sf (l+p)" V (4.12) 

y y y 

and 

1 + p = (1 + p)*. (4.14) 

Here 

R = u + — ^-(1 + p)~^~ (4.15) 

y ~ 1 

and 

S = -u + + p) 1 ^-. (4.16) 

y - i 

These variables as before in the linear case play a crucial role on the numerical, as well as physi- 
cal, boundary conditions. In this nonlinear situation, the boundary condition that we derived for the 
linear case does not hold. Thus, we settle for a condition that is applicable for nonlinear problems 
found in Thompson [9]. This is conceptually very simple once the equations are cast in the Riemann 
variables. The principle here is that at z = L , for the waves not to reflect in the computational domain 
the speed of the incoming variable must be zero. Thus, applying this condition in equation (4.12) for 
the incoming variable, we obtain 

S' = + Mu - (1 + p ) 1 ^ ](1 + p) 1 ^ - ZL^+p) 1 ^ (4.17) 

Y Y Y 

This condition will be the replacing radiation condition for the nonlinear situation. 

Numerical implementation of the system of equations (4.8) yields considerable problems. Even 
though this is a hyperbolic system, it is not an integrable system. It cannot be cast in conservation or 
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weak conservation form. For these types of systems, numerical methods are not readily available in the 
literature. Thus, we had to modify the Lax - WendrofFs method described in section 6. However, the 
numerical treatment of the boundary conditions remain the same. 

Again, the same test case that was applied to the linear system was simulated here. The source 
was taken to be the same. The results are given in Figure 5. Comparing the results obtained in Figure 
4., we see an excellent agreement between the combination of the two solutions of the linearized prob- 
lems and the full nonlinear solution. Again, we have used a surface plot to show the long time behavior 
of the solution. We note in Figure 5. that there is a slight discrepancy in the results after a long time. 
This is due to the inaccuracies in the nonlinear boundary condition (4.17). It should be noted once 
again that these solutions were obtained from a fine grid (400 grid points) and to obtain the full non- 
linear problem in higher dimension can be tremendously costly. This is where the concept of solving 
two successive linear problems becomes attractive. Indeed, we followed this principle to obtain the 
three-dimensional axisymmetric problem which is discussed in the next section. 


5. Three-Dimensional Axisymmetric Problem. 

Here we consider the procedures for solving the three-dimensional problem. In particular, we 
solve the problems (P) and (P’) as described in section 1. The procedures for obtaining the radiation 
boundary conditions are similar to that which we considered in the last section on the full one- 
dimensional nonlinear situation. The difficulty here is that the problem is multi-dimensional. Thus, for 
the purpose of obtaining radiation conditions, we consider the Riemann variables in the direction where 
we wish to impose the boundary condition. Let us consider the computational domain shown in Figure 

6. Suppose on r lt which is located at r = R , we want to obtain appropriate radiation condition. We 
apply the procedure described in section 4, by considering the system of equations only in the r direc- 
tion. The first two equations that we obtain are 


where 


and 


R t + R r 


1 

— + — w - 

y 


u_ 

» 

r 


s,-s r = 


1 

w, w + 

Y 




(5.1) 

(5.2) 


S 
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Again, by using Thompson’s procedure to obtain the radiation condition at r = R , we set the speed of 
the incoming variable 5 equal to zero, which is equivalent to suppressing the r derivative. This yields 

S, = w z - —w + — (5.3) 

Y r. 

This is exactly the boundary condition given in equation (1.17a). A similar consideration is given in 
problem (P’) to obtain the boundary condition (1.25a). 

To obtain the boundary condition at z = L , we rewrite the equations with Riemann variables in 
the z direction. The first two equations of importance are 

R t + R z = y(/> - a) + ^ - - r - u r , (5.4) 

S t - S z = •!(/>- a) - ^ + f + «r. (5.5) 

Again by the principle used earlier, we obtain the boundary condition 

„ 1 w u , r f\ 

S, - ~(p - a) - — + — + u r _ (5.6) 

A similar consideration yields the boundary condition (1.25b) for problem (P*) at z = L. 

It is important to compare the difference in the nature of the wave propagation with and without 

the gravity term. The discussion we have in sections 3 and 4 gives most of the detail of the problem 

with the gravity term. Here we want to consider the equations without the gravity term. Also, we want 
to compare the nonlinear situation that corresponds to the equations with the gravity. Again, formula- 
tion of the problem plays a crucial role. However, the principle of deriving a sequence of linear prob- 

lems remains the same. The expansions that we need here to linearize the governing equations (1.1) - 
(1.3) are as follows: 

p = 1 + ea + z 2 a x + ..., (5.7a) 

p - 1 + zp + z 2 pi + ..., (5.7b) 

and 

q = eu + e 2 ui + • ■ • . (5.8) 

It is important to notice that the perturbation parameter e has a different meaning in this situa- 
tion. In the gravity case, it was a measure of the entropy production. While the same notion is issued 
from equation (1.3), the measure of acoustic pressure can drastically vary. The first resulting linear sys- 
tem from the expansion is 

u 


G t + U r + W 2 


r 


(5.9) 
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«( + —Pr = 0 , 


W* + ~P t = 0, 


Pi + Y(“r + W f ) = -Y— + /. 


(5.10) 

(5.11) 

(5.12) 


The initial conditions are 

p=a=u = w = 0 at t = 0. (5.13) 

The boundary condition at z = 0 is 

w = 0 for t > 0. (5.14) 

The radiation conditions are derived from the principles that we used for the gravity case. They are 

(5.15a) 


(w - p), = u r + — at z = L, 

Y r 


(« “ ^7P)t = w 2 + — at r = R. 
Y r 


(5.15b) 


Similarly, the problem for the nonlinear correction can be formulated. The equations and boun- 
dary conditions are listed below: 


«i u 

a ]( + u j r + w u = [(na) r + (wc), + a—], 

r r 

«1( + yPir = ~[uu r + wu z - yap,], 

Wi, + yPiz = -[uw r + wv, - yap,], 

Mi M 

+ Y(«lr + Wlz) = ~Y— + Of ~ yp(— + U r + W z ) - (up r + wp z ). 


(5.16) 

(5.17) 

(5.18) 

(5.19) 


As before, the initial conditions are zero and the radiation conditions have the same form. At z = L , 

(W] - ypO, = ~{uw r + ww z - ywpj + [u lr + y- + p{u r + u r + w z ) + y(Kp r + wp z )].(5.20) 
and at r = R , 

(«1 - ~Pl)t - ~[uu r + wu z - — p r ] + w u + — + p(— + u r + w z ) + — ( up r + wp z ). (5.21) 

Y Y t r y 

The numerical method and the treatment of boundary conditions for this case is similar to the one 
employed in the gravity case. Details are found in the next section. 
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5.1 Numerical Results and Discussions. 

Numerical simulations are made by considering a pulse source above the ground which oscillates 
sinusoidally for a finite time as before. The form is similar to the one that we used for the one dimen- 
sional simulation. It is 


-a(r - r 0 ) 2 . _ 

e 0 sin 2k t 


f(r,t) = 


0 < (S T, Ir — r 0 1 < e 


otherwise. 


(5.22) 


Here, a = 20 and T = 1 and e = .005. Moreover, r = (r,z) and r 0 = (0,1). The time step was 
chosen to be At = .01. Calculations were carried out up to 1300 steps or 13 periods. 

The results reported in this section are three-dimensional plots of pressure distribution over the 
entire computational domain at different time intervals. The results reported in Figures 7 through 11 
are the nonlinear solutions of the gravitational field equations. Each Figure shows the progress of the 
wave at the end increment of one period. The results reported in Figures 12 through 16 correspond to 
the same situation without the gravity term. As we see from these figures, the pressure field has higher 
intensity in the absence of gravity than the case that includes gravity. However, after a long time the 
intensity of the gravitational field becomes higher. This is particularly true for large distances. Figure 
17. shows the progression of the ratio of the maximum amplitude of the waves in each case. There are 
two graphs, one is along the ground z - 0 (Solid line) and the other at z = 0.5 (dashed line). Ini- 
tially, the ratio of the pressure is higher for the no gravity case. As time progresses, the ratio decreases 
and eventually the gravitational field begins to dominate. The nondimensional distances used here in 
the z and r directions are 3 and 10 respectively. Since the lengths are nondimensionalized by the scale 
height, the distance in the r direction is roughly 100 miles. Thus, in the gravitational field, the sound 
tends to propagate longer distances, as suggested by our theory. Our computations are limited for such 
short distances due to the cost in computations. The grid size was 75x250. To do more extensive calcu- 
lations, such as longer ranges, one requires a larger computational facility. The results reported here 
took roughly 32 hours of computations on a SUN Microsystem Computer (3/260 - with floating point 
accelerator). It took 15 minutes on a CRAY - XMP-48. These considerations give a rough idea of the 
large computational requirements of the problem. 

A more realistic situation that could be considered is to include temperature variation and an 
improved model of the atmosphere. A better atmospheric model validated with experiments are avail- 
able in [8], We shall include this realistic model in a future continuation of this paper. 
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6. Numerical Schemes Used in the Computations. 

Here we briefly describe the key numerical methods that we used in our computations. The cen- 
tral method is the well-known Lax-Wendroff method. In one dimensional conservative systems this 
method is commonly known. For details, we refer to [10]. For nonconservative systems, the authors 
are not aware of any extensions. Our extension that we used for the computations are presented. Also, 
the multidimensional extensions that we used here and are not available in the literature are given here. 
We begin with the one-dimensional scheme and the numerical evaluation of the Riemann variables 
used to treat the boundary conditions. 


6.1 Scheme for Linearized One-Dimensional Equation. 

In this case, our equation takes the following form : 

Ur + F z = H(u,z) 

In the interior we use the Lax-Wendroff scheme, which is 


u i+ 


i+ 1 2 = t(«/ + u /+i) - :TT"( F (u/+i) - F (“/)) + ■^■H(—(u i ' + 1 + u :/),z i) i = 0.1..JV-1 


2Az 


a. <+£i *+a i i t+ai <+ai 

u/ +A ' = u/- — HF(u , 2 ) - F(u. , 2 )) + AfH(4(u ; 2 + u. i 2 ),z ; ) i = 1..JV-1. 


2 2 


,+ 2 


r+4i 


At t , At 


Az i + f 


2 


,+ 2 *-T 


At the boundary z = L , the Riemann variables satisfy the following conditions: 

S,= GO i), 

R t +R z = W (u). 

The corresponding discretizations for these equations (see [11]), are 

S(Uat +A ‘) = S( Uj(r) + S( u^r.j) - S( uffJi*) + 2AtQ(u \ ) 

N ~2 

(l+——)R (u at + a 0 = R ( u jv) + ^ (ujv-i) - R + 2AtW (u 2 ) 

Az N~ - 

~(R (u£) - R (u^_,) - R (U^O) 

A z 


At z=0 we have u = 0 , and we use the equation satisfied by S to obtain the numerical boundary 
condition. 

It should be pointed out that in the second order nonlinear correction, the right hand side of all 
equations contain the derivatives of the first order problem. We evaluate these quantities using the fol- 
lowing formulas : 
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«/= («'“ “‘" 4 ')t7 


,+ ¥ , ,+ f ‘“Tn 1 

U, 2 = (« 2 - M Z )-^J 


u'iz = («i+l - «/) — 


f + 41 r+ 4* r+ 1 
Mil = («.. l 2 - M. i 2 ) — 
,+ 2 “2 AZ 


6.2. Scheme for One-Dimensional Nonlinear Equations. 

In this case, our equation takes the form 

u, + A(u)F 2 = H(u) 

As mentioned earlier, we modify the Lax-Wendroff scheme to obtain a numerical scheme. In the 
event that A is linear, the scheme coincides with the Lax-Wendroff scheme. 

u'i+i = -y( u 'i+i + u'i) * = 0...JV-1. 

U;+i T = u * 1+1 - ~ F («« 0 ) + Y H ( tf *+i> z i+ i) 1 = 


t+— 1 (+ At f+ Al 

U. = ^(u 1 2 + U. I 2 ) 

2 ,+ 2 -2 


f+ Al a f f+ Ai f+ Al f+ — 

. * + A i'TT 2 \ U _ / 17 / .. 2 \ _ r/ lt 2 . A /' TT . 2 


u/ +A '= u/+ A(u, 2 )-=^(F(u. ! 2 ) - F(u ! 2 )) + A/H(ui 2 ,z,) i = 1, • • ■ N-l 
Az .+ T -5 


6.3. Scheme for Linearized 3-D Equation. 

In this case the equations have the form 

u,+F z + G r = H(u,r) 

Here again, we modify the Lax-Wendroff scheme. Our modification is similar to the Richtmyers 
scheme given in [10]. Moreover, the equations are in the cylindrical coordinate system. Thus, there is a 
singularity on the axis r=0. Our scheme and treatment of the singularity are given below: 

3 1 

u'i j = ju'f j + — «•+!,/ + u'i-1./ + u'/.y+i + 1) 


u/.y A '= u/j-|j(F(u‘,. + 1 J )-F(u' / _ u )) - Aj-(G(u‘ i>;+1 ) - G(u 
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+ H(vfij,rj)At i=l,.JV-l, j = l, ■ • • M — l 


i 


u;., 2 = + u !j At ) 


uiy = - mi-ij)) 


At / m— t + At\ /■'»/ — f+Af\\ . ¥»/ — f+At x 


2 A r 


(Wj+i) - G(u/j_{)) + H(u/j i = 1, ■ • tf-1. j = l.Jtf-1 


At r=0 there is a singularity. We evaluate the singular term using L’Hopitals rule as indicated 
below: 


u 


r 


u r 


Then the scheme is modified at the origin by the following scheme: 
(predictor step) 


3 1 

u\\ 0 = -ju‘ j>0 + -^(u , i+li o+u , i _ 1) o+ 2u' u ) 


, 2 , 

« r.\0 = T"“M 

Ar 

t 0*7 / r r \ , 03 f t t v 

11 * u/. lt0 ) + — (u/ +u - «/-u) 


u/, + o Af - uYo + AfX^o 


Where X' i>0 = X , | >0 (u*i,o,u , rii0 ,u‘ I ij0 ) 
(corrector step) 


, t+At 
l i. 0 “ 



Tr'+Ar - 
u r i,0 ” 


_2_77'+Af 

Ar U ‘'° 


Uri.o = •^■0.9(u/+ 1 ,o - Uillfi) + o.l (u£tf-u£ff)) 


„ f+Af _ Tj-r+Ar , Af^f+Af 
u *,0 ~ u i\0 + — A i,0 


u v *+ Af v -t+Af 

where X^o — X^q (u/q 


— r+Ar — r+Afx 
» Ur j t Q » U z i t o ) • 


On the boundary z=L , we use the Riemann variables R and S together with the nonreflection 
condition. 


R, + R z = H i(u,.,u) 
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S, = H 2 ( u r ,u) 


The scheme for these is as follows: 

u/rj = 0.9u hj + 0.05(u^j_, + ufij+i) 

Uat- 2 ,/ = 0.9u/r_2,y + 0.05(u^_2,y_i + u/r- 2 ,/+i) 

*N-2,j = 0.9u£. A jJy + 0.05(u/r_2,y-i + u s + - & 2,j+ 1 ) 

“r + ^ = 2^r( U JV-l,y+l “ UN-l.y-l) + ^( u iV-W + l “ “N-W-l) 

0+ “)* = * (UjV,y) + R (UW-2.;) " * (UN- A 2 !;) 

+ 2Af// 1 (u r ,u) - £(Jl (u fij) - R (u' N - 2 j) - R (u ■$.%) 

r+ Ai f+AL 

S(mI^j { ) = 5(u^r t y) + 5(u^v_2,y) - S (ujvt- 2 j) + 2Atf/ 2 (u r > u ) 


The treatment of boundary conditions at r = i? boundary are similar to that of z = L above. 
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Figure captions 

1. Comparison between solutions of (P-W) and (P-Ll) at 1000 time steps. 

2. Comparison between solutions of (P-W) and (P-Ll) at 7000 time steps. 

3. Pressure time history in the computational domain for (P-Ll) (10 periods or 10,000 time steps) ( 
e~ z zp ) 

4. Pressure time history in the computational domain for nonlinear correction (10 periods or 10,000 
time steps) ( e~ z [zp + z 2 p\] ) 

5. Pressure time history in the computational domain for full nonlinear problem (10 periods or 
10,000 time steps) 

6. Three-Dimensional computational domain. 

7. Propagation of acoustic gravity waves at time t = 1 period (100 steps) 

8. Propagation of acoustic gravity waves at time t = 3 periods (300 steps) 

9. Propagation of acoustic gravity waves at time t = 5 periods (500 steps) 

10. Propagation of acoustic gravity waves at time t = 7 periods (700 steps) 

11. Propagation of acoustic gravity waves at time t = 9 periods (900 steps) 

12. Propagation of acoustic waves without gravity at time t = 1 period (100 steps) 

13. Propagation of acoustic waves without gravity at time t = 3 periods (300 steps) 

14. Propagation of acoustic waves without gravity at time t = 5 periods (500 steps) 

15. Propagation of acoustic waves without gravity at time t = 7 periods (700 steps) 

17. Propagation of acoustic waves without gravity at time t = 9 periods (900 steps) 

18. Ratio of maximum amplitude ( no gravity/with gravity) 
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Pressure time history in the computational domain for (P-Ll) 
(10 periods or 10,000 time steps) (e" z ep) 



t ^ 


Figure 3 

Pressure time history in the computational domain for nonlinear correction 
(10 periods or 10,000 time steps) (e‘ z [ep + e 2 pj]) 
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Propagation of acoustic gravity waves at time t = 1 period (100 steps) 



Figure 7 

Propagation of acoustic gravity waves at time t = 3 periods (300 steps) 



Figure 8 
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Propagation of acoustic gravity waves at time t = 5 periods (500 steps) 

source 



Figure 9 

Propagation of acoustic gravity waves at time t = 7 periods (700 steps) 



Figure 10 
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Propagation of acoustic gravity waves at time t = 9 periods (900 steps) 



Figure 11 

Propagation of acoustic waves without gravity 
at time t = 1 period (100 steps) 



Figure 12 
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Propagation of acoustic waves without gravity 
at time t = 3 periods (300 steps) 



Figure 13 

Propagation of acoustic waves without gravity 
at time t = 5 periods (500 steps) 



Figure 14 
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Propagation of acoustic waves without gravity 
at time t = 7 periods (700 steps) 



Figure 15 

Propagation of acoustic waves without gravity 
at time t = 9 periods (900 steps) 



Figure 16 






NASA 

National Aeronautics and 
Space Administration 


1 Report No 

NASA CR-4157 


Report Documentation Page 


2 Government Accession No 


3 Recipient's Catalog No 


4 Title and Subtitle 


Linear and Nonlinear Acoustic Wave Propagation 
in the Atmosphere 


5 Report Date 
June 1988 


6 Performing Organization Code 


7 Author(s) 


S. I. Hariharan and Yu Ping 


9 Performing Organization Name and Address 

University of Akron 

Department of Mathematical Sciences 

Akron, OH 44325 

12 Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665-5225 


15 Supplementary Notes 

Langley Technical Monitor: John S. Preisser 


8 Performing Organization Report No 


10 Work Unit No 
505-61-11-02 


11 Contract or Grant No 

NAG1-624 


13 Type of Report and Period Covered 
Contractor Report 


14 Sponsoring Agency Code 


16 Abstract 

This paper describes the investigation of the acoustic wave propagation theory 
and numerical implementation for the situation of an isothermal atmosphere. A 
one-dimensional model to validate an asymptotic theory and a three-dimensional 
situation to relate to a realistic situation are considered. In addition, 
nonlinear wave propagation and the numerical treatment are included. It is 
known that the gravitational effects play a crucial role in the low frequency 
acoustic wave propagation. They propagate large distances and, as such, the 
numerical treatment of those problems become difficult in terms of posing 
boundary conditions which are valid for all frequencies. Our treatment is 
discussed in detail. Open questions are posed. 


17 Key Words (Suggested by Author(sl) 

Sound propagation 
Nonlinear propagation 
Long-range sound propagation 


18 Distribution Statement 


Unclassified - Unlimited 


19 Security Classif (of this report) 
| Unclassified 


20 Security Classif (of this page) 
Unclassified 


Subject Category 71 
21 No of pages 22 Price 

38 A03 


NASA FORM 1626 OCT 86 


N AS A-Langley, 1988 


For sale by the National Technical Information Service, Springfield, Virginia 22161-2171 















End of Document 



